library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages ------------------------------------------------------------------- tidyverse 1.3.1 --
v ggplot2 3.3.5     v purrr   0.3.4
v tibble  3.1.2     v dplyr   1.0.7
v tidyr   1.1.3     v stringr 1.4.0
v readr   1.4.0     v forcats 0.5.1
-- Conflicts ---------------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
library(spatstat)
Lade nötiges Paket: spatstat.data
Lade nötiges Paket: spatstat.geom
Registered S3 method overwritten by 'spatstat.geom':
  method     from
  print.boxx cli 
spatstat.geom 2.3-0
Lade nötiges Paket: spatstat.core
Lade nötiges Paket: nlme

Attache Paket: ‘nlme’

Das folgende Objekt ist maskiert ‘package:dplyr’:

    collapse

Lade nötiges Paket: rpart
spatstat.core 2.3-0
Lade nötiges Paket: spatstat.linnet
spatstat.linnet 2.3-0

spatstat 2.2-0       (nickname: ‘That's not important right now’) 
For an introduction to spatstat, type ‘beginner’ 
library(pracma)

Attache Paket: ‘pracma’

Das folgende Objekt ist maskiert ‘package:spatstat.core’:

    rat

Die folgenden Objekte sind maskiert von ‘package:spatstat.geom’:

    integral, quad

Das folgende Objekt ist maskiert ‘package:purrr’:

    cross
library(plyr)
----------------------------------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------------------------------------------------------------

Attache Paket: ‘plyr’

Die folgenden Objekte sind maskiert von ‘package:dplyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

Das folgende Objekt ist maskiert ‘package:purrr’:

    compact
library(plotly)
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attache Paket: ‘plotly’

Die folgenden Objekte sind maskiert von ‘package:plyr’:

    arrange, mutate, rename, summarise

Das folgende Objekt ist maskiert ‘package:ggplot2’:

    last_plot

Das folgende Objekt ist maskiert ‘package:stats’:

    filter

Das folgende Objekt ist maskiert ‘package:graphics’:

    layout
library(mltools)

Attache Paket: ‘mltools’

Das folgende Objekt ist maskiert ‘package:tidyr’:

    replace_na
library(data.table)
data.table 1.14.0 using 2 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attache Paket: ‘data.table’

Das folgende Objekt ist maskiert ‘package:spatstat.geom’:

    shift

Die folgenden Objekte sind maskiert von ‘package:dplyr’:

    between, first, last

Das folgende Objekt ist maskiert ‘package:purrr’:

    transpose
library(htmlwidgets)

Function to calculate B1 and B2 of a specific omega.


calcB2 <- function(X1w, X2w, X3w) {
  norm12 = norm(X1w-X2w,"2")
  norm13 = norm(X1w-X3w, "2")
  norm23 = norm(X2w-X3w, "2")
  N = c(norm12, norm13, norm23)
  ind <- which.max(N)
  if(ind == 1){
    B2w = X3w
  }else if(ind == 2){
    B2w = X2w
  }
  else{
    B2w = X1w
  }
  return(list("B2w" = B2w, "ind" = ind))
}

Calculate ordered random variables for a specific omega.

calcOrder <- function(X1w, X2w, X3w, B2w, ind){
  norm1 = norm(X1w-B2w, "2")
  norm2 = norm(X2w-B2w, "2")
  norm3 = norm(X3w-B2w, "2")
  X3wP <- B2w
  if (ind == 1){
    if(norm2 >= norm3){
      X1wP <- X2w
      X2wP <- X3w
    }else{
      X1wP <- X3w
      X2wP <- X2w
    }
  }else if(ind == 2){
    if(norm1 >= norm3){
      X1wP <- X1w
      X2wP <- X3w
    }else{
      X1wP <- X3w
      X2wP <- X1w
    }
  }else{
    if(norm1 >= norm2){
      X1wP <- X1w
      X2wP <- X2w
    }else{
      X1wP <- X2w
      X2wP <- X1w
    }
  }
  
  return(list("X1" = X1wP, "X2" = X2wP, "X3" = X3wP))
}

Calculate F_i for a specific Omega

calcF <- function(R2, B2w, X1wP, X2wP, X3wP){
  
  # calculate F1
  # browser()
  n <- length(R2$x)
  F1ind = numeric()
  b <- Norm(X1wP - B2w)
  for(i in 1:n){
    # browser()
    x = as.numeric(R2[i,])
    a <- Norm(x - B2w)
    if(a >= b){
      F1ind <- c(F1ind,i)
    }
  }
  F1 <- R2[F1ind,]
  #R2 <- R2[-F1ind,]
  
  # calculate F2
  F2ind = numeric()
  b <- Norm(X2wP - B2w)
  # print(b)
  if(n!=0){
    for(i in 1:n){
      if(!(i %in% F1ind)){
        x = as.numeric(R2[i,])
        a <- Norm(x - B2w)
        if(a >= b){
          F2ind <- c(F2ind,i)
        }
      }
    }
  }
  F2 <- R2[F2ind,]
  #R2 <- R2[-F2ind,]
  
  #calculate F3
  #n <- length(R2$x)
  F3ind = numeric()
  b <- Norm(X3wP - B2w)
  # print(b)
  if(n != 0){
    for(i in 1:n){
      if(!(i %in% F1ind) && !(i %in% F2ind)){
        x = as.numeric(R2[i,])
        a <- Norm(x - B2w)
        if(a >= b){
          F3ind <- c(F3ind,i)
        }
      }
    }
  }
  F3 <- R2[F3ind,]
  #R2 <- R2[-F3ind,]
  
  #calculate F4 (the rest) and F0 (the empty set) 
  F4 <- R2
  
  return(list("F1ind" = F1ind, "F2ind" = F2ind, "F3ind" = F3ind))
}

Calculate Probability of one i which is part of C

calcProb <- function(sizeOmega, ind_x, ind, allF){
  count = 0
  for(i in 1:sizeOmega){
    Fi_ind <- allF[ind,i]
    # print(x)
    # print(nrow(merge(x,Fi)))
    Fi_ind = unlist(Fi_ind)
    if(ind_x %in% Fi_ind){
      count = count + 1
    }
  }
  # print(count/sizeOmega)
  return(count/sizeOmega)
}

Calculate C for one x (F_i (x) still missing)

funcC <- function(ind_x, sizeOmega, allF, empCDF){
  res = 0
  
  for(ind in 1:3){
    # print(res)
    res = res + calcProb(sizeOmega, ind_x, ind, allF)*(empCDF[ind_x]+(ind-1))
  }
  # print(res)
  return(res/4)
}
calcAllC <- function(){
  allF = replicate(sizeOmega, calcAllF())
  n <- length(R2$x)
  Cxy <- numeric(n)
  dt <- data.table(x = rnorm(100,-5,5), y = rnorm(100,-5,5))
  empCDF <- empirical_cdf(dt, ubounds=CJ(x = xR2, y = yR2))
  empCDF <- empCDF$CDF
  for(i in 1:n){
    Cxy[i] <- funcC(i, sizeOmega, allF, empCDF)
  }
  final <- R2
  final$Cxy <- Cxy
  return(final)
}

Calculate all the X_(i), F_i for all omega

calcAllF <- function(){
  # X1w = runif(2,lowerBound, upperBound)
  # X2w = runif(2,lowerBound, upperBound)
  # X3w = runif(2,lowerBound, upperBound)
  X1w = rpois(2,4)
  X2w = rpois(2,4)
  X3w = rpois(2,4)
  
  ret <- calcB2(X1w,X2w,X3w)
  B2w <- ret$B2w
  ind <- ret$ind
  
  ret <- calcOrder(X1w, X2w, X3w, B2w, ind)
  X1wP <- ret$X1
  X2wP <- ret$X2
  X3wP <- ret$X3
  
  scriptF <- calcF(R2, B2w, X1wP, X2wP, X3wP)
  
  return (scriptF)
}

Some testing

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpgYGB7cn0NCmxpYnJhcnkodGlkeXZlcnNlKQ0KbGlicmFyeShzcGF0c3RhdCkNCmxpYnJhcnkocHJhY21hKQ0KbGlicmFyeShwbHlyKQ0KbGlicmFyeShwbG90bHkpDQpsaWJyYXJ5KG1sdG9vbHMpDQpsaWJyYXJ5KGRhdGEudGFibGUpDQpsaWJyYXJ5KGh0bWx3aWRnZXRzKQ0KYGBgDQoNCg0KRnVuY3Rpb24gdG8gY2FsY3VsYXRlIEIxIGFuZCBCMiBvZiBhIHNwZWNpZmljIG9tZWdhLiANCmBgYHtyfQ0KDQpjYWxjQjIgPC0gZnVuY3Rpb24oWDF3LCBYMncsIFgzdykgew0KICBub3JtMTIgPSBub3JtKFgxdy1YMncsIjIiKQ0KICBub3JtMTMgPSBub3JtKFgxdy1YM3csICIyIikNCiAgbm9ybTIzID0gbm9ybShYMnctWDN3LCAiMiIpDQogIE4gPSBjKG5vcm0xMiwgbm9ybTEzLCBub3JtMjMpDQogIGluZCA8LSB3aGljaC5tYXgoTikNCiAgaWYoaW5kID09IDEpew0KICAgIEIydyA9IFgzdw0KICB9ZWxzZSBpZihpbmQgPT0gMil7DQogICAgQjJ3ID0gWDJ3DQogIH0NCiAgZWxzZXsNCiAgICBCMncgPSBYMXcNCiAgfQ0KICByZXR1cm4obGlzdCgiQjJ3IiA9IEIydywgImluZCIgPSBpbmQpKQ0KfQ0KDQpgYGANCg0KDQpDYWxjdWxhdGUgb3JkZXJlZCByYW5kb20gdmFyaWFibGVzIGZvciBhIHNwZWNpZmljIG9tZWdhLg0KYGBge3J9DQpjYWxjT3JkZXIgPC0gZnVuY3Rpb24oWDF3LCBYMncsIFgzdywgQjJ3LCBpbmQpew0KICBub3JtMSA9IG5vcm0oWDF3LUIydywgIjIiKQ0KICBub3JtMiA9IG5vcm0oWDJ3LUIydywgIjIiKQ0KICBub3JtMyA9IG5vcm0oWDN3LUIydywgIjIiKQ0KICBYM3dQIDwtIEIydw0KICBpZiAoaW5kID09IDEpew0KICAgIGlmKG5vcm0yID49IG5vcm0zKXsNCiAgICAgIFgxd1AgPC0gWDJ3DQogICAgICBYMndQIDwtIFgzdw0KICAgIH1lbHNlew0KICAgICAgWDF3UCA8LSBYM3cNCiAgICAgIFgyd1AgPC0gWDJ3DQogICAgfQ0KICB9ZWxzZSBpZihpbmQgPT0gMil7DQogICAgaWYobm9ybTEgPj0gbm9ybTMpew0KICAgICAgWDF3UCA8LSBYMXcNCiAgICAgIFgyd1AgPC0gWDN3DQogICAgfWVsc2V7DQogICAgICBYMXdQIDwtIFgzdw0KICAgICAgWDJ3UCA8LSBYMXcNCiAgICB9DQogIH1lbHNlew0KICAgIGlmKG5vcm0xID49IG5vcm0yKXsNCiAgICAgIFgxd1AgPC0gWDF3DQogICAgICBYMndQIDwtIFgydw0KICAgIH1lbHNlew0KICAgICAgWDF3UCA8LSBYMncNCiAgICAgIFgyd1AgPC0gWDF3DQogICAgfQ0KICB9DQogIA0KICByZXR1cm4obGlzdCgiWDEiID0gWDF3UCwgIlgyIiA9IFgyd1AsICJYMyIgPSBYM3dQKSkNCn0NCg0KYGBgDQoNCkNhbGN1bGF0ZSBGX2kgZm9yIGEgc3BlY2lmaWMgT21lZ2ENCmBgYHtyfQ0KY2FsY0YgPC0gZnVuY3Rpb24oUjIsIEIydywgWDF3UCwgWDJ3UCwgWDN3UCl7DQogIA0KICAjIGNhbGN1bGF0ZSBGMQ0KICAjIGJyb3dzZXIoKQ0KICBuIDwtIGxlbmd0aChSMiR4KQ0KICBGMWluZCA9IG51bWVyaWMoKQ0KICBiIDwtIE5vcm0oWDF3UCAtIEIydykNCiAgZm9yKGkgaW4gMTpuKXsNCiAgICAjIGJyb3dzZXIoKQ0KICAgIHggPSBhcy5udW1lcmljKFIyW2ksXSkNCiAgICBhIDwtIE5vcm0oeCAtIEIydykNCiAgICBpZihhID49IGIpew0KICAgICAgRjFpbmQgPC0gYyhGMWluZCxpKQ0KICAgIH0NCiAgfQ0KICBGMSA8LSBSMltGMWluZCxdDQogICNSMiA8LSBSMlstRjFpbmQsXQ0KICANCiAgIyBjYWxjdWxhdGUgRjINCiAgRjJpbmQgPSBudW1lcmljKCkNCiAgYiA8LSBOb3JtKFgyd1AgLSBCMncpDQogICMgcHJpbnQoYikNCiAgaWYobiE9MCl7DQogICAgZm9yKGkgaW4gMTpuKXsNCiAgICAgIGlmKCEoaSAlaW4lIEYxaW5kKSl7DQogICAgICAgIHggPSBhcy5udW1lcmljKFIyW2ksXSkNCiAgICAgICAgYSA8LSBOb3JtKHggLSBCMncpDQogICAgICAgIGlmKGEgPj0gYil7DQogICAgICAgICAgRjJpbmQgPC0gYyhGMmluZCxpKQ0KICAgICAgICB9DQogICAgICB9DQogICAgfQ0KICB9DQogIEYyIDwtIFIyW0YyaW5kLF0NCiAgI1IyIDwtIFIyWy1GMmluZCxdDQogIA0KICAjY2FsY3VsYXRlIEYzDQogICNuIDwtIGxlbmd0aChSMiR4KQ0KICBGM2luZCA9IG51bWVyaWMoKQ0KICBiIDwtIE5vcm0oWDN3UCAtIEIydykNCiAgIyBwcmludChiKQ0KICBpZihuICE9IDApew0KICAgIGZvcihpIGluIDE6bil7DQogICAgICBpZighKGkgJWluJSBGMWluZCkgJiYgIShpICVpbiUgRjJpbmQpKXsNCiAgICAgICAgeCA9IGFzLm51bWVyaWMoUjJbaSxdKQ0KICAgICAgICBhIDwtIE5vcm0oeCAtIEIydykNCiAgICAgICAgaWYoYSA+PSBiKXsNCiAgICAgICAgICBGM2luZCA8LSBjKEYzaW5kLGkpDQogICAgICAgIH0NCiAgICAgIH0NCiAgICB9DQogIH0NCiAgRjMgPC0gUjJbRjNpbmQsXQ0KICAjUjIgPC0gUjJbLUYzaW5kLF0NCiAgDQogICNjYWxjdWxhdGUgRjQgKHRoZSByZXN0KSBhbmQgRjAgKHRoZSBlbXB0eSBzZXQpIA0KICBGNCA8LSBSMg0KICANCiAgcmV0dXJuKGxpc3QoIkYxaW5kIiA9IEYxaW5kLCAiRjJpbmQiID0gRjJpbmQsICJGM2luZCIgPSBGM2luZCkpDQp9DQpgYGANCg0KDQoNCkNhbGN1bGF0ZSBQcm9iYWJpbGl0eSBvZiBvbmUgaSB3aGljaCBpcyBwYXJ0IG9mIEMNCmBgYHtyfQ0KY2FsY1Byb2IgPC0gZnVuY3Rpb24oc2l6ZU9tZWdhLCBpbmRfeCwgaW5kLCBhbGxGKXsNCiAgY291bnQgPSAwDQogIGZvcihpIGluIDE6c2l6ZU9tZWdhKXsNCiAgICBGaV9pbmQgPC0gYWxsRltpbmQsaV0NCiAgICAjIHByaW50KHgpDQogICAgIyBwcmludChucm93KG1lcmdlKHgsRmkpKSkNCiAgICBGaV9pbmQgPSB1bmxpc3QoRmlfaW5kKQ0KICAgIGlmKGluZF94ICVpbiUgRmlfaW5kKXsNCiAgICAgIGNvdW50ID0gY291bnQgKyAxDQogICAgfQ0KICB9DQogICMgcHJpbnQoY291bnQvc2l6ZU9tZWdhKQ0KICByZXR1cm4oY291bnQvc2l6ZU9tZWdhKQ0KfQ0KYGBgDQoNCg0KQ2FsY3VsYXRlIEMgZm9yIG9uZSB4IChGX2kgKHgpIHN0aWxsIG1pc3NpbmcpDQpgYGB7cn0NCmZ1bmNDIDwtIGZ1bmN0aW9uKGluZF94LCBzaXplT21lZ2EsIGFsbEYsIGVtcENERil7DQogIHJlcyA9IDANCiAgDQogIGZvcihpbmQgaW4gMTozKXsNCiAgICAjIHByaW50KHJlcykNCiAgICByZXMgPSByZXMgKyBjYWxjUHJvYihzaXplT21lZ2EsIGluZF94LCBpbmQsIGFsbEYpKihlbXBDREZbaW5kX3hdKyhpbmQtMSkpDQogIH0NCiAgIyBwcmludChyZXMpDQogIHJldHVybihyZXMvNCkNCn0NCmBgYA0KDQpgYGB7cn0NCmNhbGNBbGxDIDwtIGZ1bmN0aW9uKCl7DQogIGFsbEYgPSByZXBsaWNhdGUoc2l6ZU9tZWdhLCBjYWxjQWxsRigpKQ0KICBuIDwtIGxlbmd0aChSMiR4KQ0KICBDeHkgPC0gbnVtZXJpYyhuKQ0KICBkdCA8LSBkYXRhLnRhYmxlKHggPSBybm9ybSgxMDAsLTUsNSksIHkgPSBybm9ybSgxMDAsLTUsNSkpDQogIGVtcENERiA8LSBlbXBpcmljYWxfY2RmKGR0LCB1Ym91bmRzPUNKKHggPSB4UjIsIHkgPSB5UjIpKQ0KICBlbXBDREYgPC0gZW1wQ0RGJENERg0KICBmb3IoaSBpbiAxOm4pew0KICAgIEN4eVtpXSA8LSBmdW5jQyhpLCBzaXplT21lZ2EsIGFsbEYsIGVtcENERikNCiAgfQ0KICBmaW5hbCA8LSBSMg0KICBmaW5hbCRDeHkgPC0gQ3h5DQogIHJldHVybihmaW5hbCkNCn0NCmBgYA0KDQpDYWxjdWxhdGUgYWxsIHRoZSBYXyhpKSwgRl9pIGZvciBhbGwgb21lZ2EgDQpgYGB7cn0NCmNhbGNBbGxGIDwtIGZ1bmN0aW9uKCl7DQogICMgWDF3ID0gcnVuaWYoMixsb3dlckJvdW5kLCB1cHBlckJvdW5kKQ0KICAjIFgydyA9IHJ1bmlmKDIsbG93ZXJCb3VuZCwgdXBwZXJCb3VuZCkNCiAgIyBYM3cgPSBydW5pZigyLGxvd2VyQm91bmQsIHVwcGVyQm91bmQpDQogIFgxdyA9IHJwb2lzKDIsNCkNCiAgWDJ3ID0gcnBvaXMoMiw0KQ0KICBYM3cgPSBycG9pcygyLDQpDQogIA0KICByZXQgPC0gY2FsY0IyKFgxdyxYMncsWDN3KQ0KICBCMncgPC0gcmV0JEIydw0KICBpbmQgPC0gcmV0JGluZA0KICANCiAgcmV0IDwtIGNhbGNPcmRlcihYMXcsIFgydywgWDN3LCBCMncsIGluZCkNCiAgWDF3UCA8LSByZXQkWDENCiAgWDJ3UCA8LSByZXQkWDINCiAgWDN3UCA8LSByZXQkWDMNCiAgDQogIHNjcmlwdEYgPC0gY2FsY0YoUjIsIEIydywgWDF3UCwgWDJ3UCwgWDN3UCkNCiAgDQogIHJldHVybiAoc2NyaXB0RikNCn0NCg0KYGBgDQoNCg0KYGBge3J9DQojc2V0LnNlZWQoNDIpDQpzaXplT21lZ2EgPC0gMjANCmxvd2VyQm91bmQgPC0gLTENCnVwcGVyQm91bmQgPC0gMQ0KDQp4UjIgPC0gc2VxKC0yMCwyMCwxKQ0KeVIyIDwtIHNlcSgtMjAsMjAsMSkNClIyIDwtIGV4cGFuZC5ncmlkKHggPSB4UjIsIHkgPSB5UjIpDQp0ZW1wIDwtIFIyJHgNClIyJHggPC0gUjIkeQ0KUjIkeSA8LSB0ZW1wDQpmaW5hbCA8LSBjYWxjQWxsQygpDQptYXRyaXhDIDwtIG1hdHJpeChmaW5hbCRDeHksIG5yb3c9IHNxcnQobGVuZ3RoKGZpbmFsJEN4eSkpLGJ5cm93PUZBTFNFKQ0KUG9pc3NvblYxIDwtIHBsb3RfbHkoeCA9IHhSMiwgeSA9IHlSMiwgeiA9IG1hdHJpeEMsIHR5cGUgPSAic3VyZmFjZSIpDQpVbmlmVjENCmpwZWcoIlBvaXNzb240LTIwXzIwX09tMjAuanBlZyIpDQpwZXJzcCh4UjIsIHlSMiwgbWF0cml4Qyxjb2wgPSAic3ByaW5nZ3JlZW4iLCBzaGFkZSA9IDAuNSwgeGxhYiA9ICJ4XzEiLCB5bGFiID0gInhfMiIsIHpsYWIgPSAiQyh4KSIpDQpkZXYub2ZmKCkNCmpwZWcoIlBvaXNzb240LTIwXzIwX09tMjBDb250b3VyLmpwZWciKQ0KY29udG91cih4UjIsIHlSMiwgbWF0cml4QykNCmRldi5vZmYoKQ0KYGBgDQoNCg0KYGBge3J9DQoNCmBgYA0KDQoNClNvbWUgdGVzdGluZw0KDQpgYGB7cn0NCg0KYGBgDQoNCg==